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Monte Carlo simulation with PENELOPE (v. 2003) is applied to calculate Leksell 
Gamma Knife® dose distributions for heterogeneous phantoms. The usual spherical water 
phantom is modified with a spherical bone shell simulating the skull and an air-filled 
cube simulating the frontal or maxillary sinuses. Different simulations of the 201 source 
(— I I configuration of the Gamma Knife have been carried out with a simplified model of the 

Qh' geometry of the source channel of the Gamma Knife recently tested for both single source and 

. multisource configurations. The dose distributions determined for heterogeneous phantoms 

^ i including the bone- and/or air-tissue interfaces show non negligible differences with respect to 

those calculated for a homogeneous one, mainly when the Gamma Knife isocenter approaches 
^ ' the separation surfaces. Our findings confirm an important underdosage (~10%) nearby the 

• i-H . air-tissue interface, in accordance with previous results obtained with PENELOPE code with 

a procedure different to ours. On the other hand, the presence of the spherical shell simulating 
' the skull produces a few percent underdosage at the isocenter wherever it is situated. 

1. Introduction 

P : 0.,.,„aPla„® ,GP, (EMta me) i. a co.p«e.base. t„at.e„t d»e p.a.,i„. sys.e„, 

designed to calculate the dose distributions of the Leksell GammaKnife® (GK) for 
stereotactic radiosurgery of certain brain diseases. As almost all radiosurgery planning 
systems, GP is quite simplistic. Using a standard set of beam data, the dose distributions in 
patients are calculated by adding those corresponding to each one of the 201 beams of the 
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CJ I GK actually present in each particular treatment (Wu et al 1990, Wu 1992). GP assumes 

c/3 ' homogeneous target media and tissue heterogeneities are not taken into account (Yu and 

Sheppard 2003). 

p |! However, in stereotactic radiosurgery, Solberg et al (1998) have pointed out a remarkable 

disagreement between Monte Carlo (MC) results and those predicted by the usual planning 
systems, in case inhomogeneous phantoms are considered. 

In the investigation of dose perturbations produced by heterogeneities, MC has showed 
5^ I up as a useful tool, mainly because it accounts, in an adequate way, for the lack of electron 

equilibrium near interfaces. For the GK, Cheung et al (2001), using the EGS4 MC code, 
have found discrepancies up to 25% in case of extreme irradiation conditions, mainly near 
tissue interfaces and dose edges. This contrasts with the sub-millimeter accuracy with which 
GK operates (Elekta 1992). 

In this paper we have investigated the effects of bone- and air-tissue interfaces on 
dosimetric calculations involving the GK. To simulate the GK, a simplified geometry model 
of the source channels is considered to perform the calculations. This model was proposed 
in Al-Dweri et al (2004) and it is based onto the characteristics shown by the beams after 
they pass through the treatment helmets. It has been shown that the collimation system of 
each source channel acts as a "mathematical collimator" in which a point source, situated 
at the center of the active core of the GK source, emits photons inside the cone defined by 
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the point source and the output helmet colhmators. If an homogeneous target phantom is 
considered, this simphfied model of the GK produces doses in agreement with those found if 
the full geometry of the source channel is considered, with those calculated by other authors 
with various MC codes and with the predictions of GP, for both a single source (Al-Dweri et 
al 2004) and different multisource (Al-Dweri and Lallena 2004) configurations. 

In this work we want to use the simplified geometry model of the GK to calculate doses in 
case of heterogeneous target phantoms, including bone- and air-tissue interfaces. Simulations 
have been performed by using the version 2003 of PENELOPE (Salvat et al 2003). We 
compare our findings for the 201 source configuration with those obtained by Cheung et 
al (2001) with EGS4 and by Moskvin et al (2004) with PENELOPE (v. 2001). Different 
situations of the GK isocenter (both far and near the interfaces) are considered. 

2. Material and Methods 

2.1. Leksell Gamma Knife® model 

To study the effect of the heterogeneities, we have used different configurations of the phantom 
depicted in figure E It is chosen to be a sphere with 80 mm of radius made of water except 
for the two shadow regions. Region i is a cube with a side of 30 mm and with its center at 
50 mm of the center of the phantom as shown in the figure. It is considered to be made of 
material mi which can be air ("a"), to simulate the maxillary or frontal sinuses, or water. 
Region 2 consists of a 5 mm width spherical shell with its external surface situated at 5 mm 
of the phantom surface. We have considered this shell made of material m2 which can be 
either bone ( "b" ) , to simulate the skull, or water ( "w" ) . The different phantom configurations 
have been labeled as Vmim2- With this notation, T^^w labels the homogeneous phantom. The 
origin of the coordinate system is chosen to be at the center of the phantom, as indicated 
in the figure. In the figures below, the different regions relevant to the calculations will be 
shown in gray scales and labeled with the corresponding number in italic. 




Figure 1. Schematic plot of the phantom used in our simulations. 
As mentioned above, each one of the 201 sources of the GK are simulated according the 
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simplified geometry model which is described in detail in Al-Dweri et al (2004). It consists of 
a point source emitting the initial photons in the cone defined by source itself and the helmet 
outer collimators. The coordinates of the 201 point sources can be found in Al-Dweri and 
Lallena (2004). They are distributed in the z < region and correspond to the situation in 
which the isocenter of the GK coincides with the center of the phantom. For the simulations 
described below, in which the isocenter is situated at different positions, these coordinates 
must be shifted. The position of the isocenter appears explicitly in the figures as l{xi,yi,zi\, 
with the values of the coordinates in mm. 

Due to the fact that the distribution of the sources is not completely uniform, no 
cylindrical symmetry is shown by the system. Thus, the doses we have calculated depend 
on the three cartesian coordinates, D*^™i™2^(x, y, z). The superscript refers to the materials 
of the particular phantom 7^11111112 considered in the simulation. Throughout the paper, the 
values of the coordinates are given in mm. 

Cheung et al (2001) used a phantom similar to our Pwb- It included the 5 mm width bone 
shell at 5 mm of the phantom surface as in our case, but with a full diameter of 180 mm. This 
phantom cannot be positioned inside the treatment helmets in such a way that the isocenter 
of the GK approaches the skull interface. That is why in our simulations we have chosen 
the phantom described above, which is slightly smaller. In any case, we have compared our 
results with those of Cheung et al by performing simulations with their phantom, which we 
label T^wb- To do that we have calculated the quantity 

which corresponds to the dose obtained for Vwh divided by its maximum. 

A first evaluation of the effects of the different interfaces has been obtained by calculating 
the relative differences 

^(wwj (^xi,yi,zi) 

between the doses obtained at the isocenter for the heterogeneous, Pmim2 > ^'iid homogeneous, 
T^wwj phantoms. 

In addition, we have calculated the quantity 

j(mim2)c y = _±L {-^^y^Z) (o) 

in order to analyze the differences observed in the dose profiles calculated for the different 
phantoms. 



2.2. Monte Carlo calculations 

PENELOPE (v. 2003) (Salvat et al 2003) has been the MC code used to perform the 
calculations. PENELOPE permits the simulation of the coupled transport of electrons and 
photons, for an energy range from a few hundred eV up to 1 GeV, for arbitrary materials. 
PENELOPE provides an accurate description of the particle transport near interfaces. 

Photons are simulated in PENELOPE in a detailed way. Electrons and positrons are 
simulated by means of a mixed scheme which includes two types of events: hard events, 
which are simulated in detail and are characterized by polar angular deflections or energy 
losses larger than certain cutoff values, and soft events, which are described in terms of a 
condensed simulation based on a multiple scattering theory (Salvat et al 2003). The tracking 
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Table 1. PENELOPE tracking parameters of the materials assumed in our 
simulations. -Eabs(7) and i?abs(e~,e+) stand for the absorption energies corresponding 
to photons and electrons and positrons, respectively. 



materials 


Air 


Bone and Water 


EMl) [keV] 


1.0 


1.0 


Sabs(e-,e+) [keV] 


0.1 


50.0 


Ci 


0.05 


0.1 


C2 


0.05 


0.05 


Wcc [keV] 


5.0 


5.0 


W„. [keV] 


1.0 


1.0 


^max [cm] 


JQ35 


2035 



Table 2. Composition of the materials assumed in the MC simulations performed 
in this work. The values correspond to the weight fraction of each element in the 
material. Also the densities are quoted. The three materials have been generated 
with the code material included in the PENELOPE package and correspond to the 
numbers 104, 119 and 277 respectively in the material database of the MC code. 





Air 


Bone 


Water 


H 




0.52790 


0.111894 


C 


0.000124 


0.19247 




N 


0.755267 


0.01603 







0.231781 


0.21311 


0.888106 


Mg 




0.00068 




P 




0.01879 




S 




0.00052 




Ar 


0.012827 






Ca 




0.03050 




density [g cm~'^] 


0.0012048 


1.85 


1.0 



is controlled by means of the four parameters Ci, C2, VFcc and VFcr, as well as the absorption 
energies. All these parameters must be fixed for the materials present in the geometry 
considered in the simulation. Tabled shows the values we have assumed in our simulations. 
In addition we have fixed the parameter Smax = 10'^^ in all the simulations performed. 

The initial source was selected by sampling uniformly between the 201 sources. Initial 
photons were emitted with the average energy 1.25 MeV and uniformly in the corresponding 
emission cone. 

The simulation geometry has been described by means of the geometrical package 
PENGEOM of PENELOPE. The three materials assumed in the simulations performed (air, 
bone and water) have been generated with the code material included in the PENELOPE 
package. Table |21 gives the composition and densities of these materials. 

To score the doses we have considered cubic voxels with Ax = Ay = Az = 1 mm, 
for the 18 and 14 mm helmets, and Az = Ay = Az = 0.5 mm, for the 8 and 4 mm ones. 
In the calculation of the doses yi, zi), used to determine the relative differences 

A^^™^^(3;i, yi, zi) as given by equation @, voxels with double width in the y direction were 
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considered. 

The number of histories followed in each simulation has been 3 • 10*^. This permitted 
to maintain the statistical uncertainties under reasonable levels. The uncertainties given 
throughout the paper correspond to la. In much of the figures, the error bars do not show 
up because they are smaller than the symbols used. 

3. Results 

3.1. Comparison with EGS4 calculations 

First of all we have compared our results with those obtained by Cheung et al (2001) using the 
EGS4 code. They considered the 18 mm helmet for two situations of the isocenter: I[0,-66,0] 
and I[0, 0,-69]. Their results (solid curves) are compared with our findings (open squares) in 
figureEl where the values of -Dnorm(^c, y, z), as given by equation (P), are plotted for the three 
cartesian axes and the two positions of the isocenter. As we can see, both calculations are in 
good agreement for the two cases considered. 

In the figure, also the predictions of the GP, quoted by Cheung et al (2001), are included 
(dashed curves). For the isocenter at I[0,-66,0], a discrepancy between MC simulations and 
GP results is observed in the far negative y region (see medium left panel). This is due to 
the fact that the GP does not take into account the interfaces and assumes all tissue to be 
uniformly represented by water. The same situation is not observed when the isocenter is at 
I[0, 0,-69] (see right lower panel), because in that case the dose is roughly zero before reaching 
the interface (for z ~ —80 mm). 

3.2. Comparison with CammaPlan predictions 

Figure shows a comparison of different simulations performed for the phantom Pww (open 
squares) with GP predictions of Hamad and Mherat (2005) (dashed curves). In the upper 
panels the isocenter is situated at I[0,34,0] and the profiles along y axis are shown for the 18 
and 4 mm helmets. The isocenter situated at I[0,66,0] in medium panels where the profiles 
along X and z axes are shown for the 14 mm helmet. Finally, the profiles along y axis are 
shown in lower panels, for the 14 and 8 mm helmets and the isocenter situated at I[0,70,0]. 

As we can see, the simulation for the water phantom produces results in very good 
agreement with the GP predictions. Below (see section [3.4(1 . the effects of the interfaces in 
these cases will be analyzed and it will be clear that GP cannot describe these effects. 

This is evident also in right tail of the profiles shown in the lower panels. We can see 
how a clear discrepancy between the simulation and the GP appears at y = 80 mm, that is, 
at the external border of the phantom. There in, an interface air-water is considered in the 
simulation, while GP does not take into account such a situation. 

3.3. Effects of the tissue inhomogeneities on dose at the isocenter 

Now we analyze the results obtained for different positions of the isocenter of the GK, paying 
special attention to the situations in which the isocenter is close to the interfaces. 

First, we have investigated the effects of tissue inhomogeneities on the doses calculated 
at the isocenter. We have varied its position by fixing the coordinate yi at different values 
ranging from -70 mm to 70 mm and maintaining xi = zi = 0. The results obtained for the 
relative differences with respect to the homogeneous phantom, A^^^\o,yi,0), as given by 
equation 0, are shown in figured Therein, upper panels correspond to the phantom P„b; 
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Figure 2. Dose profiles at tlie isocenter, relative to their respective maxima, in 
percentage, for the 18 mm helmet. Panels show the results along x (upper), y 
(medium) and z (lower) axes. Left panels correspond to the isocenter situated at 
I[0,-66,0], while in the right panels the isocenter is at I[0, 0,-69]. Open squares are the 
results of our simulations. Full curves correspond to EGS4 results quoted Cheung et 
al (2001). Dashed curves correspond to the predictions of GP quoted by the same 
authors. The phantom "P^b has been considered in this case. 



medium panels represent the results in case of the heterogeneous phantom "Paw^ and, finally, 
in the lower panels the phantom V^b has been considered. We have plotted the results for 
the 18 mm (left panels) and 8 mm (right panels) helmets. Similar results are obtained for 
the 14 and 4 mm helmets. 

As we can see, the presence of the bone spherical shell (upper panels) produces a 
reduction of the dose at the isocenter with respect to that obtained for the homogeneous 
phantom. This reduction is observed at practically any position of the isocenter, being ~ 
3% for the two helmets. A higher perturbation in the dose is observed when the isocenter 
is situated exactly at the interface bone-water, yi = ±70 mm. In this case, the reduction in 
the dose rises to 5% for the 18 mm helmet and it is even larger for the 8 mm one. It seems 



7 




I[0,70,0] 




y [mm] y [mm] 



Figure 3. Dose profiles at the isoccnter, relative to their respective maxima, in 
percentage, for different positions of the isocenter. Open squares are the results of 
our simulations performed with the phantom 'Pww Dashed curves are GammaPlan 
predictions (Hamad and Mherat, 2005). Upper panels correspond to the isocenter 
situated at I[0,34,0] and show the profiles along y axis for the 18 and 4 mm helmets. 
Medium panels show the profiles along x (left) and z (right) axes for the 14 mm 
helmet with the isocenter situated at I[0,66,0]. Lower panels are similar to the upper 
ones but for the isocenter situated at I[0,70,0] and for the 14 and 8 mm helmets. 



evident from our results that the effect of the full skull, which we simulate here by means of 
the bone spherical shell, is not negligible at all. 

To have an idea about the origin of this effect, we have evaluated the reduction in the 
dose due to the bone inhomogeneity in a very simple case. We have considered a photon 
pencil beam coming from the source and reaching the phantom and we have calculated the 
dose at the isocenter neglecting scattering photons. In case the phantom V^\) is considered, 
the dose at the isocenter is proportional to (see e.g. Berger, 1968) 

D^^^\xi, yi, Zl) (X. ( ^ ) exp [-|UaS4 - H^S2, - HhS2 - l^vfSo] , (4) 
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Figure 4 

calculated at the isocenter for heterogeneous and homogeneous phantoms (see 
equation Q). Upper panel corresponds to the phantom T-'wb; medium panel 
represents the results in case the heterogeneous phantom Paw, and in the lower panels 
the phantom Vah has been considered. Results for the 18 mm (left panels) and 8 mm 
(right panels) helmets arc shown. 



where {nen/p)w is the mass energy absorption coefficient of water at the initial energy of 
the photons, Eq, and fia_, fi^ and /Xb are the attenuation coefficients of air, water and bone, 
respectively, at the same energy. The values Sj correspond to the length of the trajectory 
segments traveled in the region i of the phantom and thus, 

r = S4 + S3 + S2 + So (5) 

is the distance from the source to the isocenter. If we consider the phantom Pww; 

D^'"^^Xl,yi,Zl) OC exp [-/XaS4 - A*w(s3 + S2 + So)] . (6) 

We are interested in the fraction of both doses, which is given by 
DW(xi,yi,Zi) r . ^ 1 

— — = exp - /^b - /^w S2 • 7 

For the photon energy considered in our simulations, Eq = 1.25 MeV, the attenuation 
coefficients can be calculated easily (see Hubbell and Seltzer, 2004) and one obtains ny, = 
0.11174 cm~^ and ^„ = 0.06323 cm~^. On the other hand the length S2 can vary from source 
to source, depending on the position of the isocenter. If the phantom is centered with respect 
to the helmet, that is if the isocenter is at I[0,0,0], S2 = 5 mm for all the sources. In this 
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case the dose ratio is 0.976, and a reduction of 2.4% is found. This is the minimum reduction 
found for all positions of the isocenter. By varying them in the interval (—70 mm, 70 mm) 
in the three directions, we sample the full volume of the phantom, S2 ranges between 0.5 cm 
and 2.7 cm and the reduction due to the bone inhomogeneity varies between 2.4% and 12.2%. 
These results indicate that, as we have found in our simulations, a few percent reduction in 
the dose at the isocenter is expected due to the bone inhomogeneity, independently of the 
position of the isocenter. 

The air-tissue interface (central panels) produces a slight increase (1-2% at most) in the 
dose at the isocenter, in comparison with that found for the homogeneous phantom, when 
it is situated far from the separation surface. When the interface is approached the relative 
difference A^^™^^ increases, the dose at the isocenter for the inhomogeneous phantom is 
~5% larger than that obtained for the homogeneous phantom and this occurs until a point 
very close to separation surface is reached. Once the isocenter is situated at this position, 
the dose calculated for the heterogeneous phantom reduces strongly with respect to that of 
the homogeneous one. This reduction is ~ 15% in the inner side and ~ 10% in the outer side 
of the air cube. 

A similar situation is observed when both interfaces are present (lower panels) except 
for a general decrease of the dose at the isocenter whatever yi is. This shift toward smaller 
doses is due, as said before, to the presence of the skull, as we can see clearly in the region 
yi < 0, far from the air interface, which shows a behavior rather similar to that plotted in 
the upper panels. One should point out the fact that the overdosage observed in the region 
around yi = 30 mm for the T-'aw phantom is largely reduced when, in addition, the bone shell 
is considered. 

1.10 I — I — , 1 1 1 — I — , 1 1 — I 




1 2 3 4 5 

distance isocenter/interface [mm] 

Figure 5. Dose ratio at the isocenter for the doses calculated with the "Paw and the 
homogeneous 'Pww phantoms, as a function of the distance of the isocenter to the 
air-tissue interface and for the 8 mm helmet. Black squares are our results. Open 
squares and circles are those of Moskvin et al (2004) for two different areas of their 
phantom. 

The results obtained near the air-tissue interfaces can be compared with those found 
by Moskvin et al (2004). In figure |SJ we show the dose ratio at the isocenter for the doses 
calculated with the T^aw and the homogeneous phantoms, as a function of the distance 
of the isocenter to the air-tissue interface. Results correspond to the 8 mm helmet. Black 
squares are our results. Open squares and circles are those of Moskvin et al (2004) for two 
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different areas of their phantom. As we can see, agreement between the different calculations 
is rather reasonable. 



3.4- Effects of the tissue inhomog entities on dose profiles 

The larger effects observed appear when the isocenter is situated nearby an air-tissue 
heterogeneity. In order to analyze in detail the dose in this situation, we have calculated 
the quantities d^'^^\x,y,z), as given by equation Q, for two positions of the isocenter: 
I[0,34,0] and I[0,66,0]. In these two positions the isocenter is at 1 mm distance from the inner 
and outer sides of the air cube (region 1), respectively. Some results are plotted in figures IHl 
and [71 



I[0,34,0] 




5 75 25 35 45 55 65 5 15 25 35 45 55 65 

y [mm] y [mm] 



Figure 6. Values of dww ™^ (0, 0) in percentage, as given by equation when 
the isocenter is situated at I[0,34,0]. The profiles in the y-axis are shown for the four 
helmets. The squares correspond to the homogeneous phantom Vww Solid curves 
have been obtained with the phantom Vwb- Dashed curves represent the results in 
case of the heterogeneous phantom 'Paw Dotted curves take into account the phantom 



Figure El shows the profiles along the y-axis for the four helmets and for the isocenter 
at I[0,34,0]. Therein the squares correspond to the homogeneous phantom "Pwwi while 
solid, dashed and dotted curves have been obtained with the phantoms Pwb) "Paw and V^h, 
respectively. If only the bone is considered (solid curves), a reduction in the plateau region 
including the maximum dose is observed. This is the same reduction previously discussed for 
the dose at the isocenter. 

On the contrary, the presence of the air-tissue interface (dashed curves) produces a 
strong reduction of the dose on the "air" side (the right side in this case) of the interface 
and an enhancement of the dose profile on the "water" side (the left side in this case) of the 
separation surface. This effects are better seen for the 18 helmet. The main effect of the 
simultaneous consideration of both interfaces (dotted curves) is to cancel the overdosage on 



11 



the left side of the interface. The results here obtained are very similar to those plotted in 
the figure 6 of the work of Moskvin et al (2004). These large differences in the dose produced 
by the air-tissue heterogeneities cannot be neglected. 




100 



Figure 7. Values of d!^™^\x , y , z) in percentage, as given by equation 0, when 
the isocenter is situated at I[0,66,0]. The profiles in the three cartesian axes are shown 
for the 14 mm helmet. The lines have the same meaning as in figure El 



Figure [7| depicts the results obtained for the 14 mm helmet when the isocenter is situated 
at I[0,66,0]. The different curves correspond to the same phantoms as in the previous figure. 
Here the profiles along the three cartesian axes are plotted. Two facts deserve a comment. 
First, it is appreciable in the upper and lower panels the reduction of the dose produced by 
the presence of the air-tissue interface (dashed and dotted curves). Also, the comparison 
between both curves gives us an idea of the additional diminution produced by the bone 
shell. Second, it is again evident (see medium panel) the strong overdosage produced if the 
air-tissue interface is not taken into account, but it is also remarkable the reduction in the 
dose observed in the region 2 when the bone shell is considered (solid curve). 

To complete our analysis, we show in figure IS] results similar to those in figure IHl but for 
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the isocenter situated at I[0,70,0], that is exactly at the bone-tissue interface. Apart from 
the reductions observed in the dose nearby the air-tissue interface, here the effect of the bone 
sheU in region 2 is, as expected, stronger than in the previous case. 

To finish we point out that, as it is observed in figures IHl and |H1 and in the medium panel 
of figure [3 the dose in the interface air-water is, in all cases, ~50% of the maximum dose 
obtained for the homogeneous phantom. This result is into agreement with the findings of 
Moskvin ei al (2004). 

4. Conclusions 

In this work we have investigated the dosimetry of the GK in case of heterogeneous phantoms 
by considering a simplified source model for the single source channels. Calculations have been 
done by using the Monte Carlo code PENELOPE (v. 2003) for the configuration including 
201 unplugged sources and for different positions of the isocenter of the GK. 

The use of the simplified model produces results for the dose profiles at the isocenter 
which are in agreement with previous calculations done with EGS4, whereas they show 
discrepancies with the predictions of the GP, mainly at the interfaces. 

In general we can say that the presence of typical tissue inhomogeneities produces 
an underdosage with respect to the results obtained when an homogeneous phantom is 
considered. This happens for almost all the positions of the isocenter of the GK. This 
underdosage can reach values larger than 10% in the vicinity of air-tissue interfaces. The 
only exception of this conclusion occurs when the isocenter is situated at a distance of 
a few millimeters of the air-tissue separation surface, where an overdosage is produced. 
However, this overdosage is very small if, in addition to the tissue inhomogeneity, also the 
bone inhomogeneity is considered. 

We have analyzed the doses deposited in phantoms including bone and air 
inhomogeneities and we have found non-negligible discrepancies with the doses obtained 
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in case of water homogeneous phantom. In this respect it is worth to mention that an air 
inhomogeneity simulating the maxillary or frontal sinuses, give rise to large modifications of 
the dose profiles. 

We have found a reasonable agreement with previous calculations performed by Moskvin 
et al with PENELOPE code in the case of the air-tissue interfaces. 

In what refers to the bone-tissue inhomogeneity representing the skull, our results show a 
~ 3% underdosage at the isocenter, with respect to the doses calculated for the homogeneous 
phantom. This effect can be observed wherever the isocenter is situated. 

The discrepancies observed between the results obtained for heterogeneous and homo- 
geneous phantoms suggest that GP predictions must be corrected in order to take care of 
the air- and bone-tissues inhomogeneities, mainly in those cases in which the interfaces are 
present nearby the target area. 
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